* Figure3.do  12/09/23

do Figure3prep

* (1) Collapse to annual X adult status
collapse (sum) deaths pop, by(ucdtitle mcdtitle age year)
keep if mcd=="Opioid"

* (2) Dec 2023 update from CDC Wonder
* (2a) CDC Wonder shows no change for any year 1999-2020, either 0-17 or 18+

* (2b) miniscule update for 18+ year 2021: from 79028 to 79037
replace deaths = 79037 if year==2021 & age=="18+"

* (2c) proportionally large update for 0-17 year 2021: from 473 to 674
replace deaths = 674 if year==2021 & age=="0-17"

* (3) Construct variable for chart
gen deathrt = deaths*100000/pop
gen rtse = sqrt((deaths/pop)*(1-(deaths/pop))/pop)
gen upper = deathrt + 1.96*100000*rtse
gen lower = deathrt - 1.96*100000*rtse
gen tstat = deathrt/(100000*rtse)
gen deathrt25 = deathrt*25
gen upper25 = upper*25
gen lower25 = lower*25

* (4) Display chart
#delimit ;
twoway (rcap lower25 upper25 year if mcd=="Opioid" & age=="0-17", lcolor(blue))
  (scatter deathrt25 year if mcd=="Opioid" & age=="0-17", msymbol(S) mcolor(blue) msize(small))
  (line deathrt year if mcd=="Opioid" & age=="18+", lcolor(black))
  (line deathrt25 year if mcd=="Opioid" & age=="0-17", lcolor(blue)),
  title("Figure 3.  Opioid fatality rates among minors and adults")
  legend(order(2 "Ages 0-17, rescaled 25X" 3 "Ages 18+" 1 "95% confidence interval") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Deaths per 100K population")
  note("Note: The mortality rate for minors is multiplied by 25 to show on the same scale with adults."
  "Adult confidence intervals are not shown because they are less than 0.5 per 100K"
  "Population and mortality source: CDC Wonder.", size(vsmall)) ;
#delimit cr

* (5) Confirm that the confidence intervals are tight for adults, as stated in the legend.
*     Note that a confidence interval of 0.5 per 100K would have rtse = (0.5/100000)/(2*1.96) = 1.3 * 10^(-6)
tabstat tstat rtse pop deaths if mcd=="Opioid", stat(min max) by(age)
